close all; clear all; path(pathdef); clc;

% add subfunctions
addpath(genpath('subfunctions/'))
addpath(genpath('data_figure3and4/'))

% color library
color_list{1} = [67 84 147]/255;        % dark blue
color_list{2} = [110 153 201]/255;      % light blue
color_list{3} = [196 40 27]/255;        % dark red
color_list{4} = [218 134 121]/255;      % light red
color_list{5} = [40 127 70]/255;        % dark green
color_list{6} = [161 196 139]/255;      % light green
color_list{7} = [218 165 32]/255;       % gold
color_list{8} = [238 221 130]/255;      % light gold
color_list{9} = [128 0 128]/255;        % purple
color_list{10} = [204 153 255]/255;     % light purple
color_list{11} = [161 165 162]/255;     % grey

%% heating data %%

heating_data_filename = 'IR_heating_data_190917_2.xlsx';

[~,sheet_info] = xlsfinfo(heating_data_filename);

% retrieve all data set 
data = {};
for ii = 1:length(sheet_info)
    [num,txt,raw] = xlsread(heating_data_filename,ii);
    data{ii} = raw;
end

% sort out the measurement information
meas_info = {};
for ii = 1:length(sheet_info)
    meas_info{ii} = strsplit(sheet_info{ii},', ');%Split string according to ","
end

%% data processing %%

% input parameters
bInfo.Cell_type         = {'AB','P1','ABa','ABp','EMS','P2','ABal','ABar','ABpl','ABpr','MS','E','C','P3'};
bInfo.ChillerTemp       = '-2';      % chiller temperature 5/6/12 degC
bInfo.Position          = 'n';      % edge(e)/nucleus(n)/side(s)
bInfo.HeatedCell        = 'P1';     % heated cell

% get an index list for the input data
table_idx_list          = find_table_idx_for_relevant_data_sets(meas_info, bInfo.Cell_type, bInfo.ChillerTemp, bInfo.Position, bInfo.HeatedCell);

% data processing
datap                   = get_relevant_cell_division_data_sets(data, meas_info, table_idx_list);

%%%%% plotting %%%%%%
close(figure(1))
fh = figure(1); 
set(fh, 'position', [100 100 [1200 400]])

% plotting parameter
bPlot.bw                = 0.12; 
bPlot.FontSize          = 20;
bLegend                 = datap.IR_current;

% plot the table chart of cell division times
plot_bar_chart_of_cell_division_times(datap, table_idx_list, bInfo, bPlot, bLegend, color_list)

%% differential readout

bPlot.FontSize              = 12;
bPlot.LineWidth             = 1;

fh = figure(456); cla()
set(fh, 'position', [100 100 [180 215]])
power_list = [0 50 75 100 125 150]; 

% IR calibration
LD_reading                  = [0 25 50 100 150 200 250]; % mA
Power_mW_before_obj         = [0 0.92 3.85 9.92 15.9 21.9 27.8]; % mW
transmission                = 0.35;
mA_to_mW                    = fitplot({LD_reading(2:end), Power_mW_before_obj(2:end)}, 'linear', 'noplot');
% conversion factor from current to power
xdata_power                 = max(transmission * mA_to_mW(power_list),0);

inversion_data_mean = cell2mat(datap.rel_avg_AB_P1_time);
inversion_data_std  = cell2mat(datap.rel_std_AB_P1_time);

plot_bar_plots_with_errorbar(xdata_power(1:end-1), inversion_data_mean(1:end-1), inversion_data_std(1:end-1), bPlot, color_list)
axisinfo('IR power (mW)', '\tau_{P1}-\tau_{AB} (min)', [-0.8 5.3], [-2.2], bPlot.FontSize)

ax1 = gca; % current axes
ax1_pos = ax1.Position; % position of first axes
ax2 = axes('Position',ax1_pos,...
    'XAxisLocation','top',...
    'YAxisLocation','right',...
    'Color','none');

set(gca, 'XTick', [0.25    0.5833    0.9167])
set(gca, 'XTickLabel', {'1.1', '2.3', '3.4'})
set(gca, 'YTick', [])
axisinfo('IR power (mW)', '', [0 1], [], bPlot.FontSize)
set(0,'defaultTextFontName', 'Helvetica')

%%  absolute division timing all

idx = 1;

avg_time = datap.avg_time{idx}(1:end);
std_time = datap.std_time{idx}(1:end);

time_offset = zeros(1,length(avg_time));

% ABa/ABp
time_offset(3) = avg_time(1);
time_offset(4) = time_offset(3);

% EMS/P2
time_offset(5) = avg_time(2);
time_offset(6) = time_offset(5);

% ABal/ABar
time_offset(7) = time_offset(3) + avg_time(3);
time_offset(8) = time_offset(7);

% ABpl/ABpr
time_offset(9) = time_offset(4) + avg_time(4);
time_offset(10) = time_offset(9);

% MS/E
time_offset(11) = time_offset(5) + avg_time(5);
time_offset(12) = time_offset(11);

% C/P3
time_offset(13) = time_offset(6) + avg_time(6);
time_offset(14) = time_offset(13);

fh = figure(555); cla()
set(fh, 'position', [300 500 [360 200]])

Cell_type_str = {};
cell_order_idx = [9 10 4 1 3 7 8 12 11 5 2 6 13 14];
for ii = 1:length(avg_time)

   newidx = cell_order_idx(ii);

   xdata = [time_offset(newidx)+[0 avg_time(newidx)] fliplr(time_offset(newidx)+[0 avg_time(newidx)])];
   ydata = [[ii ii]-0.4 [ii ii]+0.4];
   
   xdata = xdata - avg_time(2);
   
   fill(xdata, ydata, color_list{2}, 'EdgeColor', color_list{1}, 'FaceAlpha', 1)
   plot([max(xdata)-std_time max(xdata)+std_time], [ii ii], 'color', color_list{1}, 'linewidth', 1)
   
   hold on
   
   Cell_type_str{ii} = bInfo.Cell_type{newidx};
   
end

set(gca, 'YTick', [1:length(avg_time)])
set(gca, 'YTickLabel', Cell_type_str)
axisinfo('Cell cycle time (min)', '', [-50 150], [0.2 14.8], 12)

%% later stages 
idx1 = 1;   % 0 mA
idx2 = 5;   % 100 mA

avg_time_ref = datap.avg_time{idx1}(1:end);
std_time_ref = datap.std_time{idx1}(1:end);

avg_time_inv = datap.avg_time{idx2}(1:end);
std_time_inv = datap.std_time{idx2}(1:end);

fractional_avg = 100*(avg_time_inv-avg_time_ref)./avg_time_ref;
fractional_std = 100*sqrt(std_time_ref.^2 + std_time_inv.^2)./avg_time_ref;

fh = figure(10); cla();
set(fh, 'position', [300 500 [180 215]])

% AB / P1
barinfo = bar(fractional_avg(1:2));
barinfo.XData = barinfo.XData;
barinfo.FaceColor = color_list{4};
barinfo.BarWidth = 0.65;
hold on

% later cells
barinfo = bar(fractional_avg(3:6));
barinfo.XData = barinfo.XData + 2;
barinfo.FaceColor = color_list{2};
barinfo.BarWidth = 0.65;
hold on

set(gca, 'XTick', 1:6)
set(gca, 'XTickLabel', bInfo.Cell_type(1:6))
xtickangle(45)

for ii = 1:length(fractional_avg)
    if ii <= 2
        plot([ii ii], [fractional_avg(ii)-fractional_std(ii) fractional_avg(ii)+fractional_std(ii)], 'color', color_list{3}, 'linewidth', 1)
    else
        plot([ii ii], [fractional_avg(ii)-fractional_std(ii) fractional_avg(ii)+fractional_std(ii)], 'color', color_list{1}, 'linewidth', 1)
    end
end
axisinfo('', 'Fractional change (%)', [0.5 6.5], [-60 20], 12)


%% show bar plots (two cell data only)

bPlot.FontSize              = 12;
bPlot.LineWidth             = 1;

% select data
bardata                     = get_relevant_cell_division_data_sets(data, meas_info, find_table_idx_for_relevant_data_sets(meas_info, bInfo.Cell_type, bInfo.ChillerTemp, bInfo.Position, bInfo.HeatedCell));
[xdata, ydata, ydataerr]    = get_two_cell_data_sets(bardata);

% IR calibration
LD_reading                  = [0 25 50 100 150 200 250]; % mA
Power_mW_before_obj         = [0 0.92 3.85 9.92 15.9 21.9 27.8]; % mW
transmission                = 0.35;
mA_to_mW                    = fitplot({LD_reading(2:end), Power_mW_before_obj(2:end)}, 'linear', 'noplot');

% conversion factor from current to power
xdata_power                 = max(transmission * mA_to_mW(xdata),0);

fh = figure(100); cla();
set(fh, 'position', [100 100 [200 220]])
plot_bar_plots_with_errorbar(xdata_power, ydata, ydataerr, bPlot, color_list);
set(gca, 'XTick', [0 3 6])
set(gca, 'YTick', [15 25 35])
axisinfo('IR power (mW)', 'Division time (min)', [], [], bPlot.FontSize)

set(fh, 'position', [100 100 [70 110]])
set(gca, 'XTick', [4.32 4.7])
set(gca, 'YTick', [23 25 27])
set(gca, 'XTickLabel', {'AB', 'P1'})
axisinfo('', 'Division time (min)', [4.1 4.9], [23 27], 9)

%% average temperature

% base temperature (degC)
T_base          = 12.3; 

% P1 or AB?
bPlotCell       = 'P1';

%%% calibration %%%
time_AB         = @(T) 1e-15*5.425*exp(1e3*10.43./(T+273)) + 1e12*15.81*exp(-1e3*8.58./(T+273));
time_P1         = @(T) 1e-15*0.7469*exp(1e3*11.03./(T+273)) + 1e12*18.38*exp(-1e3*8.529./(T+273));

%%% conversion of IR current to temperature %%%
IR_list_mA      = linspace(0, 200, 2000);
IR_list_fit     = max(transmission * mA_to_mW(IR_list_mA),0);
IR_list_data    = max(transmission * mA_to_mW(xdata),0);

dT_at_100mA     = 9.6; 
%
if strcmp(bInfo.HeatedCell,'P1')
    % average temperature model
    AB_IR_to_T_avg  = max(mA_to_mW(IR_list_mA),0)*dT_at_100mA*(4.0/10.8)/mA_to_mW(100); 
    P1_IR_to_T_avg  = max(mA_to_mW(IR_list_mA),0)*dT_at_100mA*(6.2/10.8)/mA_to_mW(100);
    % nuclear temperature model
    AB_IR_to_T_nuc  = max(mA_to_mW(IR_list_mA),0)*dT_at_100mA*(4.1/10.8)/mA_to_mW(100);
    P1_IR_to_T_nuc  = max(mA_to_mW(IR_list_mA),0)*dT_at_100mA/mA_to_mW(100);
elseif strcmp(bInfo.HeatedCell,'AB')
    % average temperature model
    AB_IR_to_T_avg  = max(mA_to_mW(IR_list_mA),0)*dT_at_100mA*(5.7/10.8)/mA_to_mW(100);
    P1_IR_to_T_avg  = max(mA_to_mW(IR_list_mA),0)*dT_at_100mA*(4.2/10.8)/mA_to_mW(100);
    % nuclear temperature model
    AB_IR_to_T_nuc  = max(mA_to_mW(IR_list_mA),0)*dT_at_100mA/mA_to_mW(100);
    P1_IR_to_T_nuc  = max(mA_to_mW(IR_list_mA),0)*dT_at_100mA*(4.1/10.8)/mA_to_mW(100);
else
end

bPlotOption.MarkerSize  = 5;
bPlotOption.CapSize     = 0;
bPlotOption.LineWidth   = 1;
bPlotOption.FontSize    = 12;
bPlotOption.MarkerList = {'o', 'o'};

% plot comparison results
close(figure(1000));
fh = figure(1000);
set(fh, 'position', [100 100 [110 150]])
if strcmp(bInfo.HeatedCell,'P1')
    if strcmp(bPlotCell, 'P1')
        % error bands
        errorband_x = [IR_list_fit' fliplr(IR_list_fit')];
        errorband_y = [time_P1(T_base + 1.1*P1_IR_to_T_avg)' fliplr((time_P1(T_base + 0.9*P1_IR_to_T_avg)'))];
        fill(errorband_x(1:end-1), errorband_y(1:end-1), color_list{4}, 'EdgeColor', 'none', 'FaceAlpha', 0.3)
        hold on
        
        errorband_y = [time_P1(T_base + 1.1*P1_IR_to_T_nuc)' fliplr((time_P1(T_base + 0.9*P1_IR_to_T_nuc)'))];
        fill(errorband_x(1:end-1), errorband_y(1:end-1), color_list{4}, 'EdgeColor', 'none', 'FaceAlpha', 0.3)
        hold on
        
        plot(IR_list_fit, time_P1(T_base + P1_IR_to_T_avg), 'color', color_list{4}, 'linewidth', 2)
        hold on
        plot(IR_list_fit, time_P1(T_base + P1_IR_to_T_nuc), ':','color', color_list{4}, 'linewidth', 2)
        
        plotcustom(IR_list_data, ydata(:,2), ydataerr(:,2), bPlotOption, 1, 3)
        box on
    elseif strcmp(bPlotCell, 'AB')
        % error bands
        errorband_x = [IR_list_fit' fliplr(IR_list_fit')];
        errorband_y = [time_AB(T_base + 1.1*AB_IR_to_T_nuc)' fliplr((time_AB(T_base + 0.9*AB_IR_to_T_nuc)'))];
        fill(errorband_x(1:end-1), errorband_y(1:end-1), color_list{2}, 'EdgeColor', 'none', 'FaceAlpha', 0.2)
        hold on
        
        errorband_y = [time_AB(T_base + 1.1*AB_IR_to_T_avg)' fliplr((time_AB(T_base + 0.9*AB_IR_to_T_avg)'))];
        fill(errorband_x(1:end-1), errorband_y(1:end-1), color_list{2}, 'EdgeColor', 'none', 'FaceAlpha', 0.2)
        hold on
                
        plot(IR_list_fit, time_AB(T_base + AB_IR_to_T_avg), 'color', color_list{2}, 'linewidth', 2)
        hold on
        plot(IR_list_fit, time_AB(T_base + AB_IR_to_T_nuc), ':', 'color', color_list{2}, 'linewidth', 2)
        
        plotcustom(IR_list_data, ydata(:,1), ydataerr(:,1), bPlotOption, 1, 1)
        box on
    else
    end
elseif strcmp(bInfo.HeatedCell,'AB')
    if strcmp(bPlotCell, 'P1')
        % error bands
        errorband_x = [IR_list_fit' fliplr(IR_list_fit')];
        errorband_y = [time_P1(T_base + 1.1*P1_IR_to_T_avg)' fliplr((time_P1(T_base + 0.9*P1_IR_to_T_avg)'))];
        fill(errorband_x(1:end-1), errorband_y(1:end-1), color_list{2}, 'EdgeColor', 'none', 'FaceAlpha', 0.3)
        hold on
        
        errorband_y = [time_P1(T_base + 1.1*P1_IR_to_T_nuc)' fliplr((time_P1(T_base + 0.9*P1_IR_to_T_nuc)'))];
        fill(errorband_x(1:end-1), errorband_y(1:end-1), color_list{2}, 'EdgeColor', 'none', 'FaceAlpha', 0.3)
        hold on
        
        plot(IR_list_fit, time_P1(T_base + P1_IR_to_T_avg), 'color', color_list{2}, 'linewidth', 2)
        hold on
        plot(IR_list_fit, time_P1(T_base + P1_IR_to_T_nuc), ':', 'color', color_list{2}, 'linewidth', 2)
        plotcustom(IR_list_data, ydata(:,2), ydataerr(:,2), bPlotOption, 1, 1)
    elseif strcmp(bPlotCell, 'AB')
        % error bands
        errorband_x = [IR_list_fit' fliplr(IR_list_fit')];
        errorband_y = [time_AB(T_base + 1.1*AB_IR_to_T_nuc)' fliplr((time_AB(T_base + 0.9*AB_IR_to_T_nuc)'))];
        fill(errorband_x(1:end-1), errorband_y(1:end-1), color_list{4}, 'EdgeColor', 'none', 'FaceAlpha', 0.2)
        hold on
        
        errorband_y = [time_AB(T_base + 1.1*AB_IR_to_T_avg)' fliplr((time_AB(T_base + 0.9*AB_IR_to_T_avg)'))];
        fill(errorband_x(1:end-1), errorband_y(1:end-1), color_list{4}, 'EdgeColor', 'none', 'FaceAlpha', 0.2)
        hold on
        
        plot(IR_list_fit, time_AB(T_base + AB_IR_to_T_avg), 'color', color_list{4}, 'linewidth', 2)
        hold on
        plot(IR_list_fit, time_AB(T_base + AB_IR_to_T_nuc), ':', 'color', color_list{4}, 'linewidth', 2)
        plotcustom(IR_list_data, ydata(:,1), ydataerr(:,1), bPlotOption, 1, 3)
    else
    end
end
axisinfo('IR power (mW)', 'Time (min)', [-0.2 6], [15 50], bPlotOption.FontSize)
set(gca, 'XTick', [0 2.5 5])
set(gca, 'YTick', [15 25 35 45])
   

%% average temperature difference

close(figure(2222));
fh = figure(2222); cla()
set(fh, 'position', [100 100 [180 180]])

if strcmp(bInfo.HeatedCell,'P1')
    % error bands
    errorband_x = [IR_list_fit' fliplr(IR_list_fit')];
    errorband_y = [0.9*(P1_IR_to_T_avg - AB_IR_to_T_avg)' 1.1*fliplr((P1_IR_to_T_avg - AB_IR_to_T_avg)')];
    fill(errorband_x(1:end-1), errorband_y(1:end-1), color_list{2}, 'EdgeColor', 'none', 'FaceAlpha', 0.2)
    hold on
    
    errorband_y = [0.9*(P1_IR_to_T_nuc - AB_IR_to_T_nuc)' 1.1*fliplr((P1_IR_to_T_nuc - AB_IR_to_T_nuc)')];
    fill(errorband_x(1:end-1), errorband_y(1:end-1), color_list{2}, 'EdgeColor', 'none', 'FaceAlpha', 0.2)
    
    plot(IR_list_fit, P1_IR_to_T_avg - AB_IR_to_T_avg, 'color', color_list{1}, 'linewidth', 2)
    hold on
    plot(IR_list_fit, P1_IR_to_T_nuc - AB_IR_to_T_nuc, 'color', color_list{1}, 'linestyle', '--', 'linewidth', 2)
    set(gca, 'XTick', [0 3 6])
    axisinfo('IR power (mW)', '\DeltaT_{P1}-\DeltaT_{AB} (K)', [0 7], [], bPlotOption.FontSize)    
elseif strcmp(bInfo.HeatedCell,'AB')
    % error bands
    errorband_x = [IR_list_fit' fliplr(IR_list_fit')];
    errorband_y = [0.9*(AB_IR_to_T_avg - P1_IR_to_T_avg)' 1.1*fliplr((AB_IR_to_T_avg - P1_IR_to_T_avg)')];
    fill(errorband_x(1:end-1), errorband_y(1:end-1), color_list{2}, 'EdgeColor', 'none', 'FaceAlpha', 0.2)
    hold on
    
    errorband_y = [0.9*(AB_IR_to_T_nuc - P1_IR_to_T_nuc)' 1.1*fliplr((AB_IR_to_T_nuc - P1_IR_to_T_nuc)')];
    fill(errorband_x(1:end-1), errorband_y(1:end-1), color_list{2}, 'EdgeColor', 'none', 'FaceAlpha', 0.2)
    
    plot(IR_list_fit, AB_IR_to_T_avg - P1_IR_to_T_avg, 'color', color_list{1}, 'linewidth', 2)
    hold on
    plot(IR_list_fit, AB_IR_to_T_nuc - P1_IR_to_T_nuc, 'color', color_list{1}, 'linestyle', '--', 'linewidth', 2)
    set(gca, 'XTick', [0 3 6])
    axisinfo('IR power (mW)', '\DeltaT_{AB}-\DeltaT_{P1} (K)', [0 7], [], bPlotOption.FontSize)
end
